We exclude models whose AIC and/or BIC value is above a certain threshold. Due to the large number of models, many AIC/BIC values are similar. Therefore, the analysis in this file should be seen as a preliminary step for model selection which is based on both the AIC/BIC values and the independence properties of the residuals (analysed in the next Rmd file).
The output is in various different files and needs to be joined in one tibble.
# All files that were generated on ukko
vec_files = list.files(paste0(params$PATH, params$JOBID))
vec_files = vec_files[grepl("arrayjob", vec_files)]
SCRIPT_PARAMS = readRDS(paste0(params$PATH, params$JOBID, "/", vec_files[1]))[[1]]$results_list$script_params
DIM_OUT = SCRIPT_PARAMS$DIM_OUT
if (file.exists(paste0(params$PATH, params$JOBID, "/",
"d_modelselection_aic_bic.rds"))){
tt_full = readRDS(paste0(params$PATH, params$JOBID, "/",
"d_modelselection_aic_bic.rds"))
} else {
tibble_list = vector("list", length(vec_files))
for (ix_file in seq_along(vec_files)){
file_this = readRDS(paste0(params$PATH, params$JOBID, "/",
vec_files[ix_file]))
SCRIPT_PARAMS_this = file_this[[1]]$results_list$script_params
N_NOISE_PARAMS_SGT = DIM_OUT^2 + DIM_OUT*3
IX_ARRAY_JOB_this = SCRIPT_PARAMS_this$IX_ARRAY_JOB
N_MODS_this = with(SCRIPT_PARAMS_this, N_MODS_PER_CORE * N_CORES)
tibble_list[[ix_file]] =
enframe(file_this) %>%
rename(nr = name) %>%
mutate(nr = nr + (IX_ARRAY_JOB_this-1)*N_MODS_this) %>%
unnest_wider(value) %>%
unnest_wider(results_list) %>%
select(nr, params_deep_final, value_final, input_integerparams) %>%
mutate(n_params = map_int(params_deep_final, length)) %>%
mutate(n_params_sys = n_params - N_NOISE_PARAMS_SGT) %>% # SGT has (too) many noise parameters, so I want to check what happens when one uses only the number of system parameters as punishment
unnest_wider(input_integerparams) %>%
mutate(punish_aic = n_params * 2/SCRIPT_PARAMS_this$N_OBS) %>%
mutate(punish_bic = n_params * log(SCRIPT_PARAMS_this$N_OBS)/SCRIPT_PARAMS_this$N_OBS) %>%
mutate(value_aic = value_final + punish_aic) %>%
mutate(value_bic = value_final + punish_bic) %>%
select(-value_final, -starts_with("punish"))
}
tt_full = reduce(tibble_list, bind_rows)
tt_full = tt_full %>%
mutate(p_plus_q = p + q) %>%
mutate(n_unstable = kappa * DIM_OUT + k) %>%
select(nr, p_plus_q, p ,q, n_unstable, contains("value"), n_params_sys)
saveRDS(tt_full,
file = paste0(params$PATH, params$JOBID, "/",
"d_modelselection_aic_bic.rds"))
}
DIM_OUT = SCRIPT_PARAMS$DIM_OUT
tt_full %>%
arrange(value_aic) %>%
mutate_if(is.numeric, round, digits = 4) %>%
DT::datatable(extensions = c('Buttons', 'Scroller'),
options = list(dom = 'Bflirtp', buttons = I('colvis'),
deferRender = TRUE, scrollY = 350, scroller = TRUE),
filter = list(position = 'top', clear = TRUE))
## This version of bslib is designed to work with rmarkdown version 2.7 or higher.
Group by AR and MA order to check if there is a trend regarding the optimal number of unstable MA roots
tt_full %>%
select(-contains("bic")) %>%
arrange(desc(p_plus_q), desc(p), desc(q), value_aic) %>%
group_by(p_plus_q) %>%
slice(1, .preserve = TRUE) %>%
arrange(value_aic)
tt_full %>%
ggplot(aes(x = p_plus_q, y = value_aic, color = as.factor(q))) +
geom_point() -> pp
plotly::ggplotly(pp)
tt_full %>%
ggplot(aes(x = n_unstable, y = value_aic, color = as.factor(p_plus_q))) +
geom_point() -> pp
plotly::ggplotly(pp)
tt_full %>%
filter(q <= 4, p <= 4) %>%
ggplot(aes(x = p_plus_q, y = value_aic, shape = as.factor(p), color = as.factor(q))) +
geom_point() +
scale_shape_discrete(name = "AR and MA \norder (p,q)") +
scale_color_discrete(name = "") -> pp
plotly::ggplotly(pp)
tt_full %>%
filter(q <= 4, p <= 4) %>%
ggplot(aes(x = p_plus_q, y = value_aic, color = as.factor(q), shape = as.factor(p))) +
geom_point() +
scale_color_discrete(name = "MA and MA \norder (q,p)") +
scale_shape_discrete(name = "") -> pp
plotly::ggplotly(pp)
# myfunction <- function(var, string) {
# print(var)
# print(string)
# result <- paste(as.character(string),'_new', sep="")
# return(result)
# }
#
#
# labels = tibble(p = 0:4, q = 0:4)
tt_full %>%
filter(p <= 3, q <= 3) %>%
ggplot(aes(x = n_unstable, y = value_aic)) +
geom_point() + facet_grid(p~q, labeller = label_both, scales = "free_y") +
ggtitle("Dependence of AIC Value w.r.t. Number of MA Roots Inside The Unit Circle") -> pp
plotly::ggplotly(pp)
tt_full %>%
filter(p <= 3, q <= 3) %>%
ggplot(aes(x = n_unstable, y = value_aic)) +
geom_point() + facet_grid(p~q, labeller = label_both) +
labs(title = "Dependence of AIC Value w.r.t. Integer-Valued Parameters") +
theme_bw() +
theme(plot.title = element_text(vjust = 1)) +
ylab("AIC Value") + xlab("Number of MA Roots Inside the Unit Circle")-> pp
ggsave(filename = "../local_data/paper_outputs/AICdependence.eps",
plot = pp,
device = "eps")
## Saving 7 x 5 in image
plotly::ggplotly(pp)
tt_full %>%
ggplot(aes(x = p_plus_q, y = value_aic, color = as.factor(q))) +
geom_point() -> pp
plotly::ggplotly(pp)
tt_full %>%
filter(p<=3, q<=3) %>%
ggplot(aes(x = n_unstable, y = value_bic, color = as.factor(p_plus_q), shape = as.factor(p_plus_q))) +
geom_point() +
theme_bw() +
ggtitle("BIC Value w.r.t. Number of MA Roots Inside the Unit Circle") +
xlab("Number of MA Roots Inside the Unit Circle") + ylab("BIC Value") +
labs(shape = "p plus q") + labs(color = "p plus q") +
theme(legend.title = element_text(hjust = 0 ,vjust = 0)) -> pp
plotly::ggplotly(pp)
tt_full %>%
filter(q <= 4, p <= 4) %>%
ggplot(aes(x = p_plus_q, y = value_bic, shape = as.factor(p), color = as.factor(q))) +
geom_point() +
scale_shape_discrete(name = "AR and MA \norder (p,q)") +
scale_color_discrete(name = "") -> pp
plotly::ggplotly(pp)
tt_full %>%
filter(q <= 4, p <= 4) %>%
ggplot(aes(x = p_plus_q, y = value_bic, color = as.factor(q), shape = as.factor(p))) +
geom_point() +
scale_color_discrete(name = "MA and MA \norder (q,p)") +
scale_shape_discrete(name = "") -> pp
plotly::ggplotly(pp)
tt_full %>%
filter(p <= 3, q <= 3) %>%
ggplot(aes(x = n_unstable, y = value_bic)) +
geom_point() + facet_grid(p~q, labeller = label_both) +
labs(title = "Dependence of BIC Value w.r.t. Integer-Valued Parameters") +
theme_bw() +
theme(plot.title = element_text(vjust = 1)) +
ylab("BIC Value") + xlab("Number of MA Roots Inside the Unit Circle")-> pp
plotly::ggplotly(pp)
ggsave(filename = "../local_data/paper_outputs/BICdependence_bq.eps",
plot = pp,
device = "eps")
## Saving 7 x 5 in image
tt_full %>%
filter(p <= 3, q <= 3) %>%
ggplot(aes(x = n_unstable, y = value_aic)) +
geom_point() + facet_grid(rows = vars(p), cols = vars(q)) -> pp
plotly::ggplotly(pp)